function [moments,nameMoments] = collectMoments(Ey,Eyy,autoEyy,autoLagsIdx)
% Selecting the moments for estimation 

ny           = size(Ey,1);
numLags      = size(autoLagsIdx,2);
numMom       = ny + ny*(ny+1)/2 + ny*numLags;
moments      = zeros(numMom,1);
idxMom       = 0;
moments(idxMom+1:ny,1)= Ey;
idxMom       = idxMom + ny;
for i=1:ny
    for j=i:ny
        idxMom = idxMom + 1;
        moments(idxMom,1) = Eyy(i,j);
    end
end
for i=1:numLags
    if ndims(autoEyy) == 3
        moments(idxMom+1:idxMom+ny,1) = diag(autoEyy(:,:,autoLagsIdx(1,i)));
    else
        moments(idxMom+1:idxMom+ny,1) = autoEyy(:,i);
    end
    idxMom = idxMom + ny;
end

if nargout > 1
    % We construct a cell with general labels of the moments
    for i=1:ny
        nameMoments{i,1} = ['E[y_t(',num2str(i),')]'];
    end
    idxMom = ny;
    for i=1:ny
        for j=i:ny
            idxMom = idxMom + 1;
            %moments(idxMom,1) = Eyy(i,j);
            nameMoments{idxMom,1} = ['E[y_t(',num2str(i),')y_t(',num2str(j),')]'];
        end
    end
    for k=1:numLags
        for i=1:ny
            idxMom = idxMom + 1;
            %moments(idxMom+1:idxMom+ny,1) = diag(autoEyy(:,:,autoLagsIdx(1,i)));
            nameMoments{idxMom,1} = ['E[y_t-',num2str(autoLagsIdx(1,k)),'(',num2str(i),')^2]'];
        end
    end
end


end

